function [C, P]=fuzzy_knnnew(d,Cp,K,beta)
% d : distance matrix
% Cp: class labels
% beta: is the parameter for the fuzzibess in membership
% K: k nearest neightbour
%computer fuzzy nearest prototype algorithm, the C class lables should be unique and in decressed order
warning off;
[N_data N_proto]=size(d);

% Find all class labels
ClassIndex=unique(Cp); %there is a bug at here, Cp should be in decreased order
N_class=length(ClassIndex); % number of different classes  


if K==1,   % sort distances only if K>1
  % 1NN
  % Select the closest prototype
  %[tmp,proto_index]=min(d,[],2); 
  %C=Cp(proto_index);
  %1 NN assign Fuzzy membership
  tempd=d';
  [tmp, proto_index]=sort(tempd);
  C=Cp(proto_index(1,:));
  
  old_proto_index=proto_index;
  proto_index=proto_index(1:K,:);
  knn_class=Cp(proto_index);
  all_class=Cp(old_proto_index);
  knn_dist=tmp;
  %%assign initial memberships for all data point based on sorting of K closest prototypes 
  %calcualte the membership
    sum_c_dist=sum(knn_dist.^(-2/(beta-1)),1);
    membership=[knn_dist.^(-2/(beta-1))./repmat(sum_c_dist,N_class,1)];
  %find di=0 and replace membership as 1 all other center's membership=0 
    [zeroD_i zeroD_j]=find(knn_dist==0);
    membership(:,zeroD_j)=0;
    lnZero=length(zeroD_i);
    for ki=1:lnZero
      membership(zeroD_i(ki),zeroD_j(ki))=1;
    end
  initial_membership= membership;
  
else 
  % Select up to K closest prototypes
  tempd=d';
  [tmp, proto_index]=sort(tempd);
  old_proto_index=proto_index;
  proto_index=proto_index(1:K,:);
  knn_class=Cp(proto_index);
  all_class=Cp(old_proto_index);

  %for j = 1:N_data, knn_dist(:,j) = tempd(old_proto_index(:,j),j); end   %K-NN's distance
  knn_dist=tmp;

  %assign memberships in all classes
  for i=1:N_class 
       classcounter(:,:,i)=cumsum(knn_class==ClassIndex(i));
  end
 
  %%assign initial memberships for all data point based on sorting of K closest prototypes 
  %calcualte the membership
    sum_c_dist=sum(knn_dist.^(-2/(beta-1)),1);
    membership=[knn_dist.^(-2/(beta-1))./repmat(sum_c_dist,N_class,1)];
  %find di=0 and replace membership as 1 all other center's membership=0 
    [zeroD_i zeroD_j]=find(knn_dist==0);
    membership(:,zeroD_j)=0;
    lnZero=length(zeroD_i);
    for ki=1:lnZero
      membership(zeroD_i(ki),zeroD_j(ki))=1;
    end
  initial_membership= membership;

  %% Vote between classes of K neighbors 
  [winner,vote_index]=max(classcounter,[],3);
  
  %%% Handle ties
  % Set index to classes that got as much votes as winner
  equal_to_winner=(repmat(winner,[1 1 N_class])==classcounter);
 
  % set index to ties
  [tie_indexi,tie_indexj]=find(sum(equal_to_winner,3)>1); % drop the winner from counter 
  
  % Go through tie cases and reset vote_index based on their membership of each class 
  for i=1:length(tie_indexi),
   tie_class_index=find(squeeze(equal_to_winner(tie_indexi(i),tie_indexj(i),:)));
   data_index=tie_indexj(i);
   [tmpmax tmpIdx]=max(initial_membership(tie_class_index,data_index));
   fortuna=tie_class_index(tmpIdx);
   vote_index(tie_indexi(i),tie_indexj(i))=fortuna;
 end
 C=ClassIndex(vote_index)';
end      %end if

%make out put membership
for j = 1:N_data
    reorder_membership(old_proto_index(:,j),j) = initial_membership(:,j);
end %sort the order of membership as original class order.
newmembership=repmat(reorder_membership',[1 1 K]);
if K==1
 P=zeros(N_data,N_class);
  if nargout>1,
    for i=1:N_data,
      P(i,ClassIndex==C(i))=1;
    end
  end
  P=P.*newmembership;
else
  P=shiftdim(classcounter,1)./repmat(shiftdim(1:K,-1), [N_data N_class 1]);
  P=P./repmat(max(P),N_data,1);
  P(isnan(P))=0;
  P=P.*newmembership;
end
  
    
